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Abstract 

We study stochastic sandpile models with a height restriction in one and 
two dimensions. A site can topple if it has a height of two, as in Manna's 
model, but, in contrast to previously studied sandpiles, here the height (or 
number of particles per site), cannot exceed two. This yields a considerable 
simplification over the unrestricted case, in which the number of states per 
site is unbounded. Two toppling rules are considered: in one, the particles 
are redistributed independently, while the other involves some cooperativity. 
We study the fixed-energy system (no input or loss of particles) using cluster 
approximations and extensive simulations, and find that it exhibits a continu- 
ous phase transition to an absorbing state at a critical value £ c of the particle 
density. The critical exponents agree with those of the unrestricted Manna 
sandpile. 
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I. INTRODUCTION 



Sandpile models are the prime example of self-organized criticality (SOC) a control 
mechanism that forces a system with an absorbing-state phase transition to its critical point 
@^|], leading to scale-invariance in the apparent absence of parameters ||. SOC in a slowly- 
driven sandpile corresponds to an absorbing-state phase transition in a model having the 
same local dynamics, but a fixed number of particles [3,fFHl0ll . The latter class of models 
have come to be called fixed-energy sandpiles (FES). While most studies of sandpiles have 
probed the driven case P,|14|5 there is great interest in understanding the scaling properties 
of FES models as well P|TT|-|i"5| . In this paper we study FES with a height restriction. 

From the theoretical standpoint, an inconvenient feature of sandpile models is that the 
number of particles per site is unbounded. This complicates attempts to derive cluster 



approximations and continuum descriptions. In Manna's stochastic sandpile [15,16], a site 



with z > 2 particles is active, i.e., can topple, sending two particles to neighboring sites. This 
suggests restricting the number of particles per site to z = 0, 1 or 2. In this work we study 
such a model, in one and two dimensions, with the goal of establishing its critical properties. 
Analyses of FES without a height restriction reveal that they exhibit a phase transition 
between an absorbing and an active state as the particle density ( is increased beyond a 
critical value [P,I7,n|; we find the same to be true of the restricted-height models. Thus 
the restricted model has nontrivial critical behavior, and represents, due to its simplicity, 
an attractive system for further theoretical analysis. Moreover, a detailed study allows us 
to address questions of universality in sandpiles, and, more generally, of absorbing-state 
phase transitions in systems with a conserved density (T^]. The balance of this paper is 
organized as follows. In Sec. II we define the models, followed by a discussion of cluster 
approximations in Sec. III. Numerical results are analyzed in Sec. IV, and in Sec. V we 
summarize our findings. 



II. MODELS 

The models are defined on a hypercubic lattice with periodic boundaries: a ring of L 
sites in one dimension, a square lattice of L x L sites in 2d. The configuration is specified by 
the number of particles z% = 0, 1, or 2 at each site i; sites with Zi = 2 are active, and have a 
toppling rate of unity. The continuous-time (i.e., sequential), Markovian dynamics consists 
of a series of toppling events at individual sites. (Maintaining the restiction z < 2 would 
be quite complicated in a simultaneous-update scheme.) When site i topples, two particles 
attempt to move to randomly chosen nearest neighbors j and j' of i. (j and j' need not be 
distinct.) Each particle transfer is accepted so long as it does not lead to a site having more 
than two particles. The next site to topple is chosen at random from a list of active sites, 
which must naturally be updated following each event. The time increment associated with 
each toppling is At = 1/Na, where Na is the number of active sites just prior to the event. 
At is the mean waiting time to the next event, if we were to choose sites blindly, instead 
of using a list. In this way, Na sites topple per unit time, consistent with each active site 
having a unit rate of toppling. 

We consider two stochastic toppling rules. In one, the two particles released when a 
site topples move independently. Any particle attempting to move to a site harboring two 
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particles is sent back to the toppling site. (Thus an attempt to send two particles from site 
j to site k, with Z}. = 1, results in z^ — 2 and Zj = 1.) We study this independent toppling 
rule in both one and two dimensions. In the other, cooperative rule, transitions that would 
transfer fewer than the maximum possible number of particles are avoided. The cooperative 
rule is studied in one dimension only. Transition probabilities for the two rules are listed in 
Table I. 



III. CLUSTER APPROXIMATIONS 

We have derived cluster approximations for the independent toppling rule at the one-site 
(i.e., simple mean-field theory) and two-site levels. While the height restriction complicates 
the analysis of transitions, it confers the advantage of a strict limit on the the number of 
variables. (To study the unrestricted sandpile using cluster approximations one must impose 
a cutoff on the height distribution p3].) 



A. One-site approximation 

At this level of approximation there are three variables, p n , with n = 0, 1 or 2, rep- 
resenting the probability of a site having exactly n particles. It is convenient to use the 
shorthand notation p n = (n). There is only one independent variable, due to the constraints 
of normalization, (0) + (1) + (2) = 1, and of fixed density, ( = (1) + 2(2). 

We begin the analysis by enumerating, in Fig. 1, the possible transitions between states 
of a single site. Each transition requires a specific local configuration (of two or three sites, 
depending on the process), and a particular redistribution of the two particles liberated 
when the active site topples. The local configuration and the choice of redistribution are 
independent events. In the one-site approximation all joint probabilities for two or more 
sites are factorized: (ij) — > (i)(j) and (ijk) — > (i)(j)(k). 

To illustrate how transition rates are evaluated we consider some examples. The tran- 



sition — > 1 requires the initial configuration [0j |_2j, i.e., an empty site with an active 
neighbor. Exactly one of the two particles must jump to the empty site; in d dimensions 
this occurs with probability (2d— l)/2d 2 . Thus the rate of transitions 0— >1 is 

where the factor 2d represents the number of nearest neighbors. 

Consider now the transition 2 — ► 1. There are two mutually exclusive paths by which 
it can be realized. In one, both particles jump to the same neighbor (the probability for 
this event is l/4d 2 ); if the neighbor bears a single particle, then only one particle will be 
transferred, as required. Thus the initial configuration must be [2~| |~1~| and the rate for this 
path is (2)(l)/2d. In the other path, the particles jump to distinct sites (the probability 
for this is l/2<i 2 ), one of which must already have two particles, while the other must have 



fewer than two. The required initial configuration is therefore _2j |_2j 2" , where % denotes 
a site with z < 2. The rate for this path is (2 — <i _1 )(2) 2 ( / 2). Evaluating the rates for the 
remaining transitions we obtain the equations of motion: 
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> = ^±(2)[m 2 -2(0)], (1) 

|(1) = ^(2)1(0) + (2)(JZ)-(1)] (2) 

and 

l (2) = ^ (2)[2(1) " (?)2 ~ 2(2)(?)] - (3) 

After eliminating the variables (0) and (1), a simple calculation shows that the stationary- 
density of active sites is 



(2) = 2- v /5-2C (4) 
which implies ( c — 1/2 regardless of d. 



B. Two-site approximation 

The dynamical variables are now the nearest-neighbor (NN) joint probabilities (ij) with 
i,j = 0, 1, or 2. There are four independent variables, due to the symmetry (ij) = (ji) (for 
i 7^ j) and the two relations noted previously. The allowed transitions between configurations 
of a NN pair of sites are shown in Fig. 2. 



Consider, for example, the transition 00 — > 01. The initial configuration must be | 00 | [2 
its probability, in the two-site approximation, is (00)(02)/(0), where (02) /(0) represents the 
conditional probability for a NN pair in state 02, given one site in state 0. To realize the 
transition, exactly one particle must be transferred from the toppling site to its neighbor in 
the 00 pair; this occurs with probability (2d— l)/2d 2 , as before. The rate for this process is 
then given by 

(2d-l) 2 (00)(02) 
2d 2 (0) 

where the additional factor of 2d— 1 represents the number of possible locations for the 
neighbor in state 2. (Note that in the loss term for (00) this rate is multiplied by 2 to 
account for the mirror-symmetric process.) Proceeding in this manner we obtain the rates 
for each of the 17 allowed transitions. These are used to generate the equations of motion for 
the pair probabilities, which are then integrated using a fourth-order Runge-Kutta scheme. 

We find Cc = 0.75 in 1-d (just as for the unrestricted model), and ( c = 0.63 in 2-d. 
(The corresponding simulation values are 0.92965 and 0.71127, respectively, as discussed in 
the following section.) The cluster approximation predictions for the active-site density are 
compared with simulation results in Fig. 3. An interesting qualitative result of the 2-site 
approximation is that active sites are anticorrelated, i.e., (22) < (2) 2 . This is expected on 
physical grounds, since, to become active, a site must have a NN that has toppled recently. 
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C. Cooperative rule 



For the cooperative rule, the evolution equations for the probabilities (0), (1), and (2) 



are 



|(0) = -1(020) + (121), 



(5) 



d_ 
~dt 



(1) = (020) - 2(121), 



and 



(6) 



-(2) = --(020) + (121). 



(7) 



To obtain the 1-site approximation we factorize all joint probabilities. There is then only 
one independent equation, for example, 



> = "< 2 > 



^(o) 2 + (i) 2 



(8) 



In the stationary state this gives (0) = V2(l) from which it follows that 



(2) 



3 + V2 



C-(V2-i) 



(9) 



The critical density is then ( c = V? - 1 ~ 0.41421. 

The smaller value of ( c here, as compared with the independent rule, reflects the fact the 
cooperative rule tends to maximize the number of active sites generated. We show below 
that the critical density ( c of the independent model is in fact slightly lower than that of 
the cooperative one. The reason for this is not immediately apparent from the transition 
rates, but would appear to lie in subtle correlations induced by the dynamics, that are not 
evident at the 1-site level. 



IV. SIMULATION RESULTS 
A. Independent rule 

We performed extensive simulations of the height-restricted FES with independent top- 
pling rule in one and two dimensions. The initial condition is generated by distributing (L d 
particles randomly among the L d sites, avoiding occupancy of any site by more than two 
particles. This yields an initial distribution that is spatially homogeneous, and uncorrelated. 
The dynamics begins once all the particles have been placed on the lattice. The particle 
number is, of course, conserved by the dynamics. 

In one dimension we study system sizes ranging from L = 100 to 5000 sites; in two 
dimensions the system comprises L x L sites with L = 10, 20, 40,.. .,320. For each L we 
study a range of particle densities ( = N/L d . The simulations consist of N s independent 
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runs, extending to a maximum time t m . (In one dimension, for example, we used N s = 10 5 , 
t m = 4000 for L = 100, and N s = 2000, t m = 2x 10 6 , for L = 5000. In two dimensions these 
parameters varied from N s = 10 5 , t m = 1000, for L = 10, to N s = 2 x 10 4 and t m = 8 x 10 4 
for L = 320.) 

Our first task is to locate the critical density ( c ; to this end we study the active-site density 
p a (t), its second moment p„(£), and the survival probability P(t). The second moment is 
used to evaluate the ratio m(t) = P 2 a {t) / p 2 a (t) ■ Figures 4 and 5 show typical results for p a (t) 
and P(t), respectively. p a (t) relaxes to a well-defined stationary value, p a ((,L), (similarly 
for m), while the exponential decay of of P(t) allows one to extract an associated lifetime, 
r(£, L). The stationary values, p a ((,L) and m((,L) are obtained by discarding the initial, 
transient portion of the data, and performing averages over the remainder, weighted by P(t), 
which measures the effective sample size. 

In a fixed-energy sandpile of linear extent L, we can only vary £ in increments of 1/ L d . To 
circumvent this limitation, work we adopt a strategy employed in a recent study of the pair 



contact process ||20 | . Given simulation results for the stationary values of p a and m, and of 
the survival time r, for a certain system size, we form least-squares cubic fits to these data, 
permitting interpolation to arbitrary ( values within the interval studied. Thus, for each L, 
we regard p a , m, and r as a functions of a continuous variable £■ (Since the properties of a 
finite system are nonsingular, the interpolation procedure seems quite natural.) Data sets 
for m, and associated cubic fits, are shown in Fig. 6. 

A well known criterion for criticality is size-independence of order-parameter moment 



ratios, typically in the form of "crossings" of Binder's reduced fourth cumulant pl |. Moment 



ratio crossings have also proven useful for fixing the critical parameter value at absorbing- 



state phase transitions [20,22,23]. We determine the value ( cr (L,L') for which m((,L) = 
m(C, L'), for successive L values. Extrapolating these data to L — > oo yields our estimate for 
C c ; Fig. 7 illustrates the procedure. Evidently the crossing values ( cr (L,L f ) converge quite 
rapidly. In two dimensions, the crossings are well described by the form ( cr {L, L') ~ ( c +aL~ b 
where a is an amplitude and b ~ 2.72. 

Analysis of the moment-ratio crossings yields ( c = 0.92965(3) in 1-d and ( c = 0.711270(3) 
in 2-d, where the figures in parentheses denote uncertainties. For comparison, we note the 
values for the unrestricted version of the model: 0.94885(7) in 1-d, 0.71695(5) in 2-d. Thus 
the height restriction yields a rather small shift in ( c , by about 2% in one dimension, and 
0.8% in 2-d. This is reasonable since, in the unrestricted model (near its critical point), 
only a small fraction of the sites have z > 2. The critical values of the moment ratio 
are: m c = 1.1596(4) in 1-d, and 1.347(2) in 2-d. While these differ significantly from 
the corresponding values for the directed percolation (DP) universality class [1.1735(5) and 
1.3257(5) in 1-d and 2-d, respectively p2| ], the moment ratios for the two classes are very 
similar. 

In studies of absorbing-state phase transitions pi, including fixed-energy sandpiles 



PJl8|[, it is common to determine the critical point by seeking a power-law dependence 
of the order parameter (p a in the present instance) and the relaxation time on the system 
size L. The former is governed by 

Pa ((,L)=L-^n(L 1 ^A) , (10) 



as expected on the basis of finite-size scaling [p5 |. (Here A = ( — ( c , and 1Z is a scaling 
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function.) Thus at the critical point (A = 0) we expect p a (( c ,L) ~ L w 1/ - L ; for the lifetime 
one has t(( c ,L) ~ L v \\l v ^. 

With ( c in hand, we may verify the power-law dependence of the order parameter and the 
lifetime on system size, as in Eq. (PH[), by interpolating the simulation data to the critical 
value Cc- Figure 8 shows that p a indeed has a power-law dependence on L; a similar plot 
(not shown) yields the same conclusion for r. From the data for the four largest systems, 
we then obtain (via least-squares linear fits), the exponent ratios (3/v±_ and v\\/v± listed in 
Table II. (The uncertainties reflect two contributions: one due to the uncertainty of the fit, 
the other, dominant one, due to the uncertainties in the values of p a and r for each L. The 
latter includes the effects of uncertainty in ( c .) 

To determine the exponent f3 we analyze the results for p a in the portion of the supercrit- 
ical regime where the graph of In p versus In A follows a power law. In two dimensions this 
procedure yields (3 = 0.661(3), 0.661(2), 0.654(3), and 0.655(2) for L = 20, 40, 80, and 160, 
respectively, leading to an estimate of (3 = 0.656(5). Fig. 9, a scaling plot of LP/ U± p a ((, L) 
versus L 1 ^ U± A for various system sizes, shows a good data collapse, verifying the finite-size 
scaling hypothesis for the order parameter, and yielding u± = 0.85. 

In one dimension it turns out that no power laws are seen if we use ( c = 0.92965 as 
determined from the FSS analysis described above. Quite clean power law dependence of p a 
is observed, however, if we use an L-dependent effective critical point Q >L in the analysis. 
We determine C, c ,l by optimizing the linearity of In p a as a function of In A, and maximizing 
the number of data points that may reasonably be fit by the power law. (For L = 1000 for 
example, we are able to fit 15 points with a correlation coefficient of 0.99996.) The resulting 
values of Cc,l and (3 are listed in Table III. Extrapolating the effective critical densities to 
infinite L (via linear regression versus L~ 1 ^ u± ) yields ( c = 0.9298(4), consistent with our 
our estimate based on moment-ratio crossings. A similar extrapolation gives (3 = 0.412(4) 
We determine the exponent u± via a data-collapse analysis, as in the 2d case (see Fig. 10). 
We obtain a good data collapse for l/u± in the range 0.60 to 0.62, leading to the estimate 
v ± = 1.64(4). 



B. Cooperative rule 

We performed extensive simulations of the cooperative model in one dimension with 
system sizes again ranging from L = 100 to L — 5000. In this case we prevented the system 
from falling into an absorbing configuration by maintaing at least two active sites. (If there 
are only two active sites, transitions which decreases the number of active sites are not 
permitted. Actually, there is only one transition of this sort, 020 — > 101.) The density of 
active sites p is then always > 2/L. But since the stationary value of p a at the critical 
point is ~ L~@/ UJ -, with /3/u± ~ 1/4, this should have a minimal effect on critical properties. 
Figure 11 shows the stationary active-site density as a function of the particle density for 
several values of L. 

We first analyzed the stationary critical properties of the model by means of the finite 
size scaling relation, Eq. fllPp. The critical density was obtained by plotting p a versus L 



for several ( values, as shown in Fig. 12. (As before, values of p a for densities between 
those accessible for a given L were obtained via interpolation.) Using the criterion of power 
law dependence of the order parameter on system size, we find ( c = 0.9788(1) and (3/v± = 
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0.245(5). As an alternative determination of ( c we used moment-ratio crossings. Figure 13 
shows the moment ratio m as a function of ( for L = 2000 and L = 5000. The two curves 
cross at ( c = 0.9788, confirming the previous result. 

Having obtained the critical particle density, we used it to find the critical exponent (3 
governing the order parameter. Figure 14 is a log-log plot of p a versus ( — ( c for several 
values of L. The slope of the straight line fitted to the data points for L = 5000 gives 
(3 = 0.417(1). From this, and our previous result for /3/v±, we obtain v±_ = 1.70(4). 

We also performed time-dependent simulations at the critical density, to measure the 
growth of the number of active sites. Here, each trial began with just one active site. For 
a given particle density (, this was realized by placing a particle at each of £L — 1 distinct 
sites, chosen at random. One of these sites was then selected randomly, and another particle 
placed there, rendering it active. In a lattice of size L = 10000, we performed from 5000 
to 6000 trials of this kind, to determine the mean number of particles, n(t), averaged over 
all trials (including those that fall into an absorbing configuration prior to time t). At the 
critical point, and for a sufficient large system, n(t) is expected to increase asymptotically 
as a power law 

n{t) ~ f (11) 
where the exponent rj is related to the exponent z — by the scaling relation TTT1. 




(12) 

in d dimensions. Our data for n(t) at ( c = 0.9788 do in fact follow a power law, and yield 
the estimate rj = 0.330(5). Using our previous result for f3/v± we then obtain z = 1.54(5). 



V. DISCUSSION 

We studied the scaling behavior of fixed-energy sandpiles that follow a stochastic dynam- 
ics similar to that of the Manna model, but with a height restriction Z{ < 2. Both versions 
of the model (i.e., the independent and cooperative toppling rules), exhibit a continuous 
phase transition between an absorbing state and an active one at a critical particle density 
Cc- One- and two-site cluster approximations do not yield very accurate predictions for the 
critical density (as is to be expected), but they correctly predict the continuous nature of 
the transition. 

As shown in Table II, the critical exponents for the present models appear to be the 
same as for the unrestricted case. In fact, there is excellent agreement between the exponent 
values for the restricted and unrestricted models, except for the exponent z = v\\/v±_ in one 



dimension. As was noted in Ref. [|1^], however, obtaining a reliable estimate for the exponent 
z from simulations is quite difficult in one dimension. (In two dimensions the estimates for 
z are in excellent agreement. It appears that the relaxation dynamics is anomalous in one 
dimension, as suggested in 

Studies of a reaction model suggest a common universality class, distinct from that of 
directed percolation, for absorbing-state phase transitions in which the order parameter is 
coupled to a second field that relaxes diffusively in the presence of activity [O]. Manna's 
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stochastic sandpile falls in this category, with the local particle density C( x )^) playing the 
role of the second field. Our results show that height restrictions and perturbations of 
the toppling rule do not alter the critical exponents, if they preserve the above-mentioned 
features, supporting universality in critical behavior far from equilibrium. 
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TABLES 



Transition 


Probability 




Independent 


Cooperative 


UzU - 


v i ni 
-> 1U1 


1 /o 

1/2 




1 /o 
1/2 




-> 2UU 


1 1 A 
1/4 




1/4 




nno 
-> UUz 


1 //l 
1/4 




1 1 A 

1/4 


i on 
Izll - 


-» zlll 


1 /o 
1/2 




1 /o 
1/2 




, i no 
-> lUz 


1/4 




1 /o 

1/2 




-> 210 


1/4 







220 - 


-> 202 


1/4 




1 




-> 211 


1/2 









-> 220 


1/4 







121 - 


-> 202 


1/2 




1 




-> 112 


1/4 









-> 211 


1/4 







122 - 


-» 212 


3/4 




1 




-> 122 


1/4 








Table I. Transition probabilities for the independent and cooperative toppling rules in one 
dimension. The transition probabilities are symmetric under reflection. 



Model 


Cc 


P/u ± 




P 


Indep. Id 


0.92965(3) 


0.247(2) 


1.45(3) 


0.412(4) 


Coop. Id 


0.9788(1) 


0.245(5) 


1.54(5) 


0.417(1) 


Unrst. Id 


0.94885(7) 


0.239(11) 


1.66(7) 


0.42(2) 


Indep. 2d 


0.711270(3) 


0.774(3) 


1.572(7) 


0.656(5) 


Unrst. 2d 


0.71695(5) 


0.78(2) 


1.57(4) 


0.64(1) 



Table II. Critical parameters of restricted and unrestricted sandpiles. Figures in parentheses 
denote uncertainties. Results for the unrestricted models are from Refs. [0 (1-d) and |T7 
(2-d). 



L 


Cc,L 


h 


500 


0.9256(1) 


0.465(3) 


1000 


0.9273(1) 


0.441(4) 


2000 


0.92815(5) 


0.431(4) 


2000 


0.92845(5) 


0.423(4) 


CO 


0.9298(4) 


0.412(4) 



Table III. Effective size- dependent critical density and apparent exponent (3 L , and L — > oo 
extrapolated values, for the 1-d model with independent toppling rule. 
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FIGURE CAPTIONS 

FIG. 1. Transitions between states of a single site, 'x' denotes a forbidden transition; 
diagonal entries are irrelevant. 

FIG. 2. Transitions between configurations of a NN pair of sites; 'x' denotes a forbidden 
transition; diagonal entries are irrelevant. 

FIG. 3. The stationary active-site density in the two-dimensional restricted-height sandpile 
as predicted by the one-site and two-site approximations, compared with the simulation 
result. 

FIG. 4. Evolution of the active-site density in the one-dimensional sandpile with height 
restriction (independent rule). L = 1000; £ = 0.93. 

FIG. 5. Evolution of the survival probability P(t), for the same parameters as in Fig. 4. 

FIG. 6. Stationary moment ratio, m((), in the one-dimensional model (independent rule). 
Squares: L = 500; +: L = 2000. Curves are cubic fits to the data. 

FIG. 7. Moment-ratio crossing values ( cr versus reciprocal system size in 1-d (independent 
rule). 

FIG. 8. Stationary active-site density p a versus system size L at the critical point in 1-d 
(+) and 2-d (squares) (independent rule). 

FIG. 9. Scaling plot of the stationary active-site density in the 2-d (independent rule). 
System sizes L = 20 (filled squares), 40 (open squares), 80 (x) and 160 (+). 

FIG. 10. Scaling plot of the stationary active-site density in 1-d (independent rule). System 
sizes L = 500 (o), 1000 (o), 2000 (x) and 5000 (+). 

FIG. 11. Stationary active-site density in 1-d versus particle density, for various system sizes 
(cooperative rule). 

FIG. 12. Stationary active-site density versus system size in 1-d (cooperative rule). 

FIG. 13. Stationary moment ratio in the one-dimensional model (cooperative rule). 

FIG. 14. Stationary active-site density versus ( — ( c for various system sizes (cooperative 
rule, 1-d). 
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